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ABSTRACT 

in 

This document describes analytical and numerical techniques for the generation of Gaussian 
density fields, which represent cosmological density perturbations. The mathematical techniques 
involved in the generation of density harmonics in fc-space, the filtering of the density fields, and 
the normalization of the power spectrum to the measured temperature fluctuations of the Cosmic 
Microwave Background, are presented in details. These techniques are well-known amongst 
experts, but the current literature lacks a formal description. I hope that this technical report 
will prove useful to new researchers moving into this field, sparing them the task of reinventing 
the wheel. 
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Gaussian density field play a major role in cosmology. There is now strong evidence that the large- 
scale structure of the universe originated from the growth, by gravitational instability, of primordial density 
fluctuations. Observations of the temperature fluctuations of the Cosmic Microwave Background (CMB) 
; I ■ indicate that these primordial fluctuation were Gaussian. 

The formation of evolution of large-scale structure in the universe is a complex problem that requires a 
numerical approach. Typically, one creates a realization of the density fluctuations at early time, and uses 
a numerical algorithm to evolve this fluctuation all the way to the present (or to some redshift of interest). 
Since it is impossible to simulate the entire universe (which might very well be infinite) , we normally assume 
that the universe is periodic at large scales. We can then divide the universe into identical cubes of volume 
Vbox = -^box' an d we onr y need to simulate one cube. This approximation is valid as long as the box size 
Lbox is much larger than any existing large-scale structure in the universe. We can rephrase this by saying 
that the box must contain a "fair" sample of the universe. 

The density field can be represented in two different ways. In the Eulerian Representation, the box is 
divided into N x N x N cells, and the density contrast S is calculated at the center of each cell. In the 
Lagrangian Representation, N p x N p x N p equal-mass particles are laid down on a cubic grid inside the box, 
and are then displaced in order to represent the density fluctuation. The choice of representation depends 
on the particular algorithm that will be used to evolve the system from these initial conditions. 



Subject headings: cosmology: theory — methods: numerical 



INTRODUCTION 



1.1. The Power Spectrum 



A Gaussian density field can be represented as a superposition of plane waves of wavevectors k and 
complex amplitudes <5£ ont , where the superscript "cont" stands for continuous, indicating that all values of 
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k are allowed. The amplitudes are related to the power spectrum P(k) by 

|C nt | 2 ocP(fc), (1) 

where k = |k|. However, there is a lot of confusion in the literature about the constant of proportionality 
between |<5£ ont | 2 and P(k), and even the units of (5£ ont can vary from one author to another. We will clarify 
this issue in §2. 

Note: Equation (1) is a convenient simplification. As we will discuss later, in a truly Gaussian random 
field, the amplitudes 5^ are determined only in a statistical sense. Their real and imaginary parts are sepa- 
rately given by a Gaussian distribution whose variance is proportional to P(k). However, using equation (1) 
greatly simplifies the derivation presented in §2, without affecting the results. 



2. THE AMPLITUDES OF THE DENSITY MODES 

We assume that the universe is periodic over a comoving cubic volume Vb ox = L^ ox . The density 
contrast S can be decomposed into a sum of plane waves. 

6 ^ = pEf Ce " ik,r - ( 2 ) 

k 

where r is the comoving position. The wavevectors k are given by 

k = (I, to, n)k , ?, m,n = —oo, —1,0, 1, oo . (3) 

where the fundamental wavenumber is 

ko = ^. (4) 

^box 

The superscript "disc" indicates that the amplitudes 5^ lsc form a discrete spectrum, that is, they are defined 
for particular, discrete values of k. 1 The factor 1/N 3 is not necessary at this point, and could be absorbed 
into the definition of <5£ 1SC - We introduce it to make the notation consistent with §3. To make 5(r) real, the 
coefficients <5£ 1SC must satisfy the reality condition: 

S-t = ($k SC T • (5) 

Our first challenge is to relate the discrete sum in equation (2) to the continuous sum of modes present 
in the real universe, and to express the amplitudes S^ lsc in terms of the power spectrum P(k). To do so, we 
will consider the rms fluctuation of the density at a certain scale R, and match the expressions obtained in 
the discrete and continuous limits. 



2.1. The Discrete Limit 

Consider a sphere of radius R centered at r . The mass inside that sphere is given by 

M(r ) = / p[l + 8(r)} d 3 r = pV sph + f d 3 rJ2 S^e^ , (6) 

Jsph(r ) Jsph(r ) k 



1 Other values of k would not satisfy the periodic boundary conditions 
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where p is the average density, V^ph = 47ri? 3 /3 is the volume of the sphere, and the integral is computed over 
that volume. The relative mass excess in the sphere is given by 



AM M{r )-pV sph 1 



M 



pV S ph 



sph Jsph(ro) 



disc — ik-r 
k e 



(7) 



We introduce the following change of variables, 



r = r + y . 



(8) 



In y-space, the sphere is now located at the origin, and equation (7) becomes 

AM 



AM, N 1 f , 

"M~ (ro) = /vn^T / 
' i w v sp h j S ph(o) k 



disCg— ikrog— iky 



(9) 



We now square this expression, and get 



/AM 



M 



16tT 2 R 6 N 6 



sph(O) 



disc — ikr — iky 



sph(O) k , 



(10) 



The variance of the density contrast at scale R is obtained by averaging the above expression over all possible 
locations Tq of the sphere inside the computational box, 



4 



AM V 



M J V hoxJVb: 

I Vbox 



V box 16tt 2 #W« J Vbo 
The integral over Vbox reduces to 



f d *r [ <?y [ rf 3, EE ,d isc ,d ;sCe -i, ye -ik'. Ze -i(k + k0. r c . (n) 

JV box Jsph(O) ^sph(O) k k , 



/ d 3 r oe - 4 ( k+k '' ro = V box 5 k ,-k' • 

J Vbox 



(12) 



We substitute this expression in equation (11), and use the Kronecker 5 to eliminate the summation over k'. 
Equation (11) reduces to 



2 _ f \ ^ I rdisc 1 2 

° R ~ 167^6^6 1^ I 
k 



-i 2 



d 3 ye' lWy 



sph(O) 



(13) 



where we used equation (5) to get iJ^^ 130 = (5 k lsc )*(5 k lsc = |<5 k lsc | 2 . The remaining integral can be evaluated 
easily (see Appendix A). Equation (13) reduces to 



4 = ^El^ isc l 2 ^ 2 (^), 



(14) 



where 



W{u) = — (sinu — ucosw) . 



(15) 



-4 - 



2.2. The Continuous Limit 

The real universe is of course not periodic, in which case all values of k are allowed. To convert the 
expressions derived in §2.1 from the discrete limit to the continuous one, consider any function /t that is 
summed over all allowed values of k. In the discrete limit, we have 



E-fdisc \ ^ rdisc ^ \ ^ rdisc7,3 

Vbox V-^ rdisc f i3i„ 

k allV.E. ^allV.E. ^ all V.E. •/ V.E. 



(16) 



where "V.E." represent a volume element in k-space, which is a cube of volume /cq centered around an 
allowed value of k (with ko = 27r/Lbox)- Assuming that the function /£ lsc does not vary significantly over 
one volume element, we can pull it inside the integral, 

k anv. E y v - E - 

Of course, the effect of integrating over the volume element, and then summing over all volume elements, is 
to effectively integrate over all fc-space, so equation (17) reduces to 

E/k isC -7^//k diSCd3fc - ( 18 ) 

We can rewrite this expression as 

E/k isc « / fZ ont d 3 k, (19) 
k J 

where the continuous and discrete functions are related by 

rcont ^box /.disc /on"! 

/k ~(2^ /k ' (20) 
Using these formulae, we can rewrite equation (2) as 

S(r) = -^ J d 3 fcC nt e- ik r , (21) 

where 

rcont ^box cdisc /00^ 

dk ~(2^ dk ■ (22) 

Let us now convert equation (14) into an integral, as we did for equation (16). We get 

< = J^bI d3kKiSCl2W2{kR) - (23) 



We substitute equation (22) into equation (23), and get 

^ = ^E J d^Sr^ikR). (24) 
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We then need to related a\ to the power spectrum P(k). The relation is given by Bunn & White (1997) as 

1 f°° 

a 2 R = — dkk 2 P(k)W 2 (kR). (25) 
21" Jo 

This relation is obtained by performing an integration over angles, using the fact that P(k) is a function of 
k = |k| only. We can "undo" this integration, simply by dividing equation (25) by An. We get 

** = (2^ / d3kP (k)W 2 (kR)- (26) 
Comparing equations (23), (24), and (26), we get 

^) = ^l^ isc | 2 = ^6lC nt | 2 - (27) 

This gives us the relation between P(k), <5£ 1SC , and <5£ ont . Notice the P(k) is neither the square of S^ lsc nor 
the square of 5£ ont . Both P{k) and (*>£ ont have dimensions of a volume while <5£ ISC is dimcnsionless. 

Equations (25) and (26) define the normalization of the power spectrum. When using any power 
spectrum obtained from the literature, it is essential to check that these relations are satisfied. Variations 
by factors of 2n between different papers are quite common. 



3. DIRECT CALCULATION OF THE DENSITY HARMONICS 

Using the formalism described in §2, we now want to compute the density harmonics 5^. We first lay 
down inside the computational volume Vb ox a regular cubic grid of size N x N x N, with grid spacing 
A = Lb ox /N. The coordinates r of the grid points are given by 

r = (a,/3, 7 )A, a, (3, 7 = 0, 1, . . . , N - 1 . (28) 

The presence of a grid results in a discretization of space, which in turns modifies the structure of the fc-space. 
In equation (3), the values of k form an infinite cubic grid in fc-space, since the indicies I, m, n can take any 
integer value from — oo to +co. However, the discretization of space limits the number of modes. Consider 
a mode with wavenumber 

k' = k+(uN,vN,wN)k , (29) 
where u, v, and w are integers. The exponential in equation (2) becomes 

e -ik'-r _ — ik-r — i[uiV,i;JV,ii>JV]-[Q,/3,7]fcoA _ e -i\t-r e -i[uN ,vN ,wN]-[a,f3,~f](2-n: /N) 

Hence, the modes k' and k are inseparable. They represent a plane wave with the same effective wavenumber. 
Consequently, we will consider a finite fc-space, in which the values I, m, n do not run from — oo to +oo, but 
instead are limited to the range to N — 1. Hence, in fc-space, the density harmonics <5(k) are also defined 
on a regular cubic grid of size N x N x N . The wavevectors are given by 



k = (I, m, n)k , I, m, n = 0, 1, . . . , N — 1 . 



(31) 
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In doing so, we are simply neglecting high-frequency modes. This is justified since the discreteness of the 
grid in r-space prevents us from resolving any structure that these modes represent. Hence, by using a grid 
in r-space, we are effectively performing a filtering of the density fluctuation at the scale of A, and such 
filtering eliminates high-frequency modes. 

The Fourier transform and inverse Fourier transform are given respectively by 

5(k) = ]T<5(rK k - r , (32) 

r 

^ = ^E^ k ) e ~ 4k ' r - ( 33 ) 

k 

Notice that equation (33) is the same as equation (2), with a slight change of notation: <5£ 1SC — > S(k). We 
introduce the notation 

5(k) - 5 lmn , (34) 
5(r) = d a01 . (35) 

Equation (32) becomes 

N-l N-l 

Simn = ^2 5 c , i g 7 e <(27r/AJV)A[ '' m '" I ' [a ' /3 ' 71 = ^2 S al 3 1 e 27rila/N e' 27 ' iml3/N e 27Tin ' y/N 

a,/3,7= a,/3,7= 

E( 27rla . 2-7r/a\ / 27rm/3 . 27rra/?\ / 27rnj . 27rn7\ 
l cos +, sin j ( cos +zsin j +Jsm j . (36 ) 

a,/3,7=0 V 7 V 7 V 7 

After expansion, this expression becomes 

&lmn — (^eee ^eoo ~t~ &oeo ~t~ ^ooe) "i" ^eeo "i" ^eoe "i" ^oee "i" $000) j (^) 



where we define 



^eee — 



J eeo — 



J eoe — 



J oeo — 



E 1 . 2-rrla 2-Km(3 2ixn r y 

<W 7 COS COS — — COS — ^- , (38) 

a,/3,7=0 

E 1 . 27r/a 27rrn/3 . 271717 

<W 7 cos cos — ^— sin — ^- , (39) 

a,/3,7=0 
JV — 1 

E 27r/a . 2-Km(3 2imy 
<W 7 cos sin — — cos , (40) 

a, 13,1=0 
N — l 

<Wr cos sin — ^— sin , (41 ) 

a,/3, 7 =0 
N—l 

E . 2irla 2irm(3 2imy 
6 af3l sin cos — — cos , (42) 

q,/3, 7 =0 

E 1 27tZq! 27rm/3 . 27rri7 

<5 Q /3 7 sm cos — ^— sin , (43) 

a, /3,7=0 
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5 v-^ 1 r . 2irla . 27TTO/3 2irnj 

6 ooe = - b aPl sin sin — cos — — , (44) 

a,/3,7=0 
TV — 1 

s c • 27tZq! . 27rm/3 . 27rn7 

a,/3,7=0 

Consider now the mode SN-i,m,n- We replace I by TV — I in equation (36), and get 

N-l N-l 
? _ \ 1 r 2iri(N-l)a/N 2irimf)/N 2irin~f/N _ \ ~- r 27ria -2nila/N 2niml3/N 2nin-y /N (Aa\ 

°N—l,mn / "aPyV e e — / VaflyV e e e . l^UJ 

a,/3,7— a,/3,7= 

Since a is an integer, the first exponential is always unity, and equation (46) reduces to 



JN-l,mn 



cos __ _ lsm __ j ^ cos _____ +ism _____ j ^ cos __ + lsm __j . (47) 



E 



_,/3,7=0 



Comparing equations (36) and (47), we see that the effect of replacing I by N — I amounts to a change of 
sign of the first sine function. In equations (38)-(45), that sine appears only in the j's for which the first 
subscript is "o." Hence, these <5's will change sign, and equation (37) will become 

+ &eoo &oeo $ooe) ~\~ ^{^eeo ~t~ 3eoe ^oee 3ooo) , (^8) 

We can directly generalize to the other indicies, or combinations of them. Replacing m by N — m 
changes the sign of the iS's for which the second subscript is "o," and replacing n by N — n changes the sign 
of the J's for which the third subscript is "o." Hence, 

Sl,N-m,n — {Seee - Seoo + S oeo $ooe) i(&eeo $eoe ~l~ $oee &000) , (^^) 

&lm,N—n = {^eee ~ Seoo — Soeo + Sooe) + i{~ $eeo ~\~ &eoe ~\~ &oee &000) 1 (50) 

$N—l,N—m,n — {&eee ~ Seoo — Soeo + S ooe ) + i(5 eeo — S eoe — 5 oee + 6 ODO ) , (51) 

^N—l,m,N—n = {&eee ~ " Seoo ~\~ S eo ~ &ooe) ~\~ i{~ &eeo ~\~ ^eoe — ^oee ~\~ &000) 1 (52) 

Sl,N—m,N—n = {$eee ~\~ Seoo ~~ Soeo — 5 ooe ) + l( — 5 eeo — 6 eoe + S oee + 8 oo) 1 (53) 

^N—l,N—m,N—n — {^eee ~\~ Seoo ~\~ S eo ~\~ &ooe) ~\~ *( — &eeo ~~ ^eoe — ^oee ~ &000) ■ (54) 

Hence, 8 different, but related harmonics can be represented by various combinations of 8 numbers: 5 eee , 
Seeo, S eoe , Soee, 4oo, ^oeo, &ooe, and 5 000 . This implies that the Fourier transform of a real field defined on a 
grid N x N x N can be represented by N 3 real numbers stored on a similar grid, even though the Fourier 
transform <5(k) is complex. Notice also that the 8 "related" modes are located, in fc-space, at the verticies of 
a rectangular box centered on the center of the fc-space grid (I, m,n = N/2). This is illustrated in Figure 1. 

From equations (37) and (48)-(54), we see that related modes form 4 pairs of complex conjugates: 

,N—m,N—n J 

(55) 

,N-m,n i (56) 
,m,N-n i (57) 

(58) 





= °N 


3lm,N-n 


= S*N 


fil,N—m,n 


= % 


N—m,N—n 


= d N 
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Fig. 1.— Representation of the system in /c-space. The 3 axes correspond to the 3 indicies Z, to, n, which 
run from to N — 1, with TV = 64 in this case. The dots indicate the particular mode (I, to, n) = (14, 15, 8) 
and the seven related modes. Notice that only one of these modes is located in the first octant, indicated by 
the dashed lines. 

3.1. General Case 

Consider first the modes for which the indicies I, to, and n are neither nor N/2. In equations (55)-(58), 
the amplitudes <5 are provided by the power spectrum. From equation (27), we get 

1/2 

(59) 

Equation (59) provides the correct normalization of the power spectrum. However, there are two problems 
with this equation. First, it provides no mean of determining the phases of the complex numbers S^ lsc , and 
second, choosing random phases would result in a field that is not Gaussian. In a truly Gaussian field, 
equation (59) is only valid in a statistical sense. The correct approach is to compute the real and imaginary 



disc | 



= N 6 



m 

Vbox 
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parts of <5£ 1SC independently, by drawing them from a Gaussian distribution, 



Re<S£ isc 



Im^ isc 



Gi(0,l)iV 3 
G 2 (0,1)7V 3 



P(k) 



2V hox 
P(k) 



-,1/2 



"I 1/2 



2V hc 



(60) 
(61) 



where Gi(0, 1) and C?2(0, 1) are random numbers drawn from a Gaussian distribution with mean and 
standard deviation 1. This ensures that the resulting field is Gaussian. 

Equations (60) and (61) provide us with the left-hand-sides of equations (55)-(58). By equating these 
expressions with equations (37), (49), (50), and (53), and considering the real and imaginary parts separately, 
we get 8 equations, 



&eee + &eoo + &oeo + Sooe = Rl , (62) 

Seeo + &eoe + ^oee + <>ooo = h , (63) 

Seee ~ ^eoo ~ ^oeo + Sooe = R2 , (64) 

~&eeo + fieoe + ^oee ~ ^000 = ^2 > (65) 

Seee — $eoo + ^oeo — Sooe = R3 7 (66) 

&eeo &eoe ~t~ ^oee ^000 ^3 3 (67) 

&eee ~l~ ^eoo ^oeo &ooe Ri , (68) 

-Seeo — &eoe + 5oee + <5 000 = I4 , (69) 



where 



Ri = Re £; TO „ , 

h = Im (5; m „ , 

i?2 = R^ &lm,N —n 1 

1 2 = lmSl m ^N- n , 
R3 = Rc5l,N~m,n, 

1 3 = lmSl t N-m,n, 
i?4 = Rec^jV-m.iV-n , 

7 4 = Im<5;,Ar_ m ,AT_„ . 



(70) 
(71) 
(72) 
(73) 
(74) 
(75) 
(76) 
(77) 



Altogether, we get two separate systems of 4 equations and 4 unknowns. Written in matrix form, these 
systems are: 



M 



" Seee 




'Ri' 




$eeo 




' h' 


Qeoo 




R2 


M 


Seoe 




-h 


Voeo 




R 3 


$oee 




h 


_ $ooe - 




. Ri . 




_ $ooo - 




-h. 



(78) 
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where the matrix M is given by 



M 



1 1 

-1 1 

1 -1 

-1 -1 



The inverse of the matrix M is simply M 1 = M/4. Hence, the S's are given by 



M 

T 



Ri 

R 2 

i?4 



ceo 
eoe 



M 

T 



h 
-h 
h 
-h 



(79) 



(80) 



3.2. On a Face 

Consider now the case where one of the indicies, say I, is equal to either or N/2. These cases correspond 
to values of k located on a face of the first octant in fe-space. The problem becomes simpler. With I = 0, 
equation (36) reduces to 



JV-l 



2-aim0/N 2win~/ /N 



JV-l 



= E * 



q/3 7 COS 



2nmj3 
N 



+ i sin - 



2Ttm(3\ ( 2irn-f . . 2-Knj\ 



COS ■ 



a, 0,7=0 

After expansion, this expression becomes 

Somn = {See + &oo) + «(^eo + S oe ) , 



N 



+ % sin - 



-)■ 



(81) 



(82) 



where 



JV-l 



E2nmP 27rn7 
Octfj-y COS — — COS ■ 



a, 13,1=0 
JV-l 



N N ' 



Seo = 



5oe = 



E2nm(3 2im^ 
<W 7 COS Sin , 

a,/3,7=0 



JV-l 



E2irm(3 2im^ 
<W 7 sin cos -jj- , 

a, 13,1=0 

E. 2-Kmj3 . 2-nn r ) 
o Q/ 3 7 sm — 77 — sin ■ 



a,/3,7=0 



N N 



(83) 
(84) 
(85) 
(86) 



Consider now the mode So.N-m.n- We replace m by N — m in equation (81), and get 

JV-l JV-l 



6 ,N-m,n= £ 5a ^i{N- m )0,N e ^ inl ,N = £ ^ , _. , 
a,/3,7= a,/3,7= 



2nip p—2Triml3 /N e 2nin7 /N 



(87) 
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Since j3 in an integer, the first exponential is always unity, and equation (87) reduces to 

N-l , 
SoM-m,n = ^ <5 Q , ) 3 7 ( 
q,/3,7=0 



27rm/3 2itm(3\ ( 2irnj 2irn r y\ 
cos — — i sm — — — cos — — h i sin ■ 



N N \ N N 



Comparing equations (81) and (88), the only difference is a change of sign of the first sine function. In 
equations (83)-(86), that sine appears only in the <5's for which the first subscript is "o." Hence, these <5's 
will change sign, and equation (82) will becomes 

So,N-m,n = (See - S oa ) + i(8 eo - S oe ) , (89) 

Similarly, we can easily show that replacing n by N — n results in a change of sign of the j's for which the 
second subscript is "o." Hence, we get 

Som.N-n = (^ee - S ao ) + i(-5 eo + 5 oe ) , (90) 
SoM-m,N-n = {See + ^oo) + i(S eo ~ $oe) ■ (91) 

These 4 modes are located at the verticies of a rectangle in /c-space, as shown in Figure 2. They form two 
pairs of complex conjugates, 

Somn = So,N-m,N-n ! (92) 

Som,N-n = So,N-m,n ■ (93) 

By equating these expressions with equations (82) and (90), and considering the real and imaginary parts 
separately, we get 4 equations, 

See+S 00 = J2l , (94) 
See + $oe = h, (95) 
See-Soo = R*, (96) 
-S eo + S oe = h , (97) 

where 

(98) 
(99) 
(100) 
(101) 

Again, the numbers R\, I\, R2, and I2 are determined from equations (60) and (61). The solutions are 

See = (102) 
Soo = (103) 

Seo = ^y^, (104) 

Soe = ^i^-- (105) 



Ri 


= Re Somn , 


h 


= Im S 0mn , 


R2 


= Re Som,N-n 


h 


= lm5om,N~n 



1 



N-l 



Fig. 2. — Representation of the system in fc-space. The dots indicate the particular mode (I, m, n) = (0, 15, 8) 
and the three related modes. 

Consider now the case I = N/2. Equation (36) reduces to 

N-l 

0N/2,mn = ^ °<*P~t e e e 
a,/3,7=0 

E 1 r , „.„ / 27rm/3 . . 2itm[3\ ( 2irn~f . . 2-im^\ 
<Wr(-l) I cos— — + zsin — — 1 I cos— — +ism— — 1 . (106) 

a,/3,7=0 ^ / V 1 / 

After expansion, this expression becomes 

Somn = {See + &oo) + i{Ko + L) , (107) 

where 

x x ( na 27rm/3 2im~f 

See = 2^ <Wy( _1 ) COS — COS — — , (108) 

ct,/3 : 7=0 
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Fig. 3. — Representation of the system in &:-space. The dots indicate the particular mode (l,m,n) = 
(32, 15,8), where 32 = N/2, and the three related modes. 

x x ( t\a 2nmf3 . 2-imj 

5 eo = 5 a p 7 (-l) "cos — sin—, (109) 

q,/3 : 7=0 

5 oe = 2^ <Wy( _1 ) sin — cos — — , (110) 

a,/3,7=0 

jf X I 1\a • 27rm ^ • 27m T ,,111 

a,/3,7=0 

We find the same expressions as for the case I = 0, the only differences being the extra factor of (— l) a 
in the definitions of the <5's. Hence, the solutions (102)-(105) are still valid in this case. These modes are 
shown in Figure 3. 

We have only considered the cases when the first index, I, is either or N/2, but these results can be 
generalized to the other indicies m and n as well, since the entire problem has cubic symmetry. 
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3.3. On an Edge 

Consider now the case when two of the indicies, say I and ra, are equal to either or N/2. These cases 
corresponds to values of k located on an edge of the first octant in fc-space. With I = m = 0, equation (36) 
reduces to 

S 00n = W^" 7/Ar = £ ^(cos^+isin^p) . (112) 

a,/3,7=0 a,/3,7=0 ^ ' 

<$00n = <5 e + i(5 G , (113) 



This expression becomes 



where 



N-l 

S e = E <W 7 cos— — , (114) 

a,/3,7=0 

d Q = 2^ <W 7 sm— — . (115) 

a,/3,7=0 

Consider now the mode <5oo,jv-n- We replace nby N — n in equation (112), and get 

N-l N-l 

6oo,N- n = E <W 2 ™ (N ~ n)7/Ar = E W^" 2 ™ 7 ^ ■ (116) 

a,/3,7=0 a,/3,7= 



Since 7 in an integer, the first exponential is always unity, and equation (116) reduces to 

N-l , 
So,0,N-n = E < ^ Q < a T ( ' 



a,/3,7=0 

Comparing equations (112) and (117), the only difference is a change of sign of the sine function. In 
equations (114) and (115), that sine appears only in the expression for <5 Q . Hence, equation (113) becomes 

Soo,N- n = 5 e - iS , (118) 

These two modes are complex conjugates, 

<W = <5oo,jv-„ ■ (119) 

They are shown in Figure 4. By equating equation (113) and (119), and considering the real and imaginary 
parts separately, we get 

5 e = Rl, (120) 

So = h, (121) 



where 



fli = Re<5 00 „, (122) 
h = ImSoon- (123) 
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Fig. 4. — Representation of the system in fc-space. The dots indicate the particular mode (I, m, n) — (0, 0, 8) 
and its related mode. 

The numbers R\ and I\ are determined from equations (60) and (61). 

Consider now the case I — N/2, m — 0, the case I = 0, m = N/2, and the case I = m = N/2. We can 
easily show that the solutions (120) and (121) are still valid, using the same approach as in §3.2. Again, 
the only difference will be extra factors of (—1)" in the definitions of S e and <5 . These modes are shown 
in Figure 5. Using the cubic symmetry, we can then show that the solutions (120) and (121) applies to all 
cases for which two of the three indicies I, m, n are equal to or N/2 (12 combinations). 



1 



N-l 



Fig. 5. — Representation of the system in fc-space. The dots indicate the particular modes (l,m,n) = 
(0, 32, 8), (32, 0, 8), and (32, 32, 8), where 32 = N/2, and their related modes. 



Finally, we consider the cases when all indicies are equal to or N/2. These cases corresponds to values 
of k located in a corner of the first octant in fc-space. With I = m = n = 0, equation (36) reduces to 

N-l 



where the subscript u stands for "uniform," since the mode 000 corresponds to a null wavenumber, or an 
infinite wavelength. Notice that since <5 Q /3 7 is real, 5 U is real as well. Hence, for this mode, there is no 
imaginary part, and S u is determined form equation (60). This mode is shown in Figure 6. 

Consider now the case I = N/2, m = 0, n = 0. Again, we can easily show that the solution (124) is 
still valid, using the same approach as in §§3.2 and 3.3. The only difference will be an extra factor of (— 1) Q 
in the expression for S u . Using the cubic symmetry, we can then show that the solution (124) applies to all 



3.4. 



In a Corner 




U ! 



(124) 



a,/3,7=0 
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Fig. 6.— Representation of the system in &;-space. The dot indicates the particular mode (I, to, n) = (0, 0, 0). 
For that particular mode, we ignore equation (124) and set S u — instead. 

cases for which the three indicies I, m, n are equal to or N/2 (8 combinations). 

Notice that there is a fundamental difference between the mode <5ooo an d the other 7 modes, 6 o.n/2, 
So. N/2,0, ■ ■ -j $n/2,n/2,n/2, shown in Figures 6 and 7, respectively. The mode <5ooo represents a perturbation 
of infinite wavelength, that is, a constant. Clearly that constant must be zero, otherwise the mean value of 
6(r) integrated over the entire volume would be nonzero, and this would violate the assumption that the 
volume contains a fair sample of the universe. Therefore, for the mode Sooo, and that mode only, we do not 
compute the amplitude |<5ooo| from the power spectrum, but instead set that amplitude equal to zero. 



3.5. Putting it All Together 

We can now count the number of independent quantities necessary to represent all the density harmonics. 
Consider first the modes Si mn for which the indicies are neither nor N/2. Excluding these values, each 
index can take N — 2 different values, which gives us (N — 2) 3 modes. As we showed in §3.1, these modes come 
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Fig. 7. — Representation of the system in /c-space. The dots indicate the particular modes (l,m,n) located 
at corners of the first octant, excluding the mode (0,0,0). For these modes, equation (124) applies. 

in groups of 8, and within each group the complex values of the harmonics can be expressed as combinations 
of 8 real numbers 5 eee , S eeo , . . ., S OOQ . Hence, it takes a total of (N — 2) 3 real numbers to represent these 
modes. 

Consider now the modes for which one of the indicies is equal to or N/2. There are 6 possibilities, 
and for each one, the remaining two indicies can take N — 2 values each (all values except and N/2). This 
gives us 6(N — 2) 2 modes. As we showed in §3.2, these modes come in groups of four, and within each group 
the complex values of the harmonics can be expressed as combinations of four real numbers S ee , 5 eo , S oe , and 
5 00 . Hence, it takes a total of 6(N — 2) 2 real numbers to represent these modes. 

Next, consider the modes for which two of the indicies are equal to or N/2. There arc 12 possibilities, 
and for each one, the remaining index can take N — 2 values (all values except and N/2). This gives us 
12(N — 2) modes. As we showed in §3.3, the amplitude these modes form complex conjugates pairs, and each 
pair can be represented by two real numbers 5 e , and 6 . Hence, it takes a total of 12(7V — 2) real numbers 
to represent these modes. 
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Finally, consider the modes for which all three indicies are equal to or N/2. There are 8 such modes. 
As we showed in §3.4, these modes are real, hence it takes 8 real numbers to represent them. The total 
number of variables necessary to represent all the density harmonics is therefore 

(N - 2f + 6(N - 2f + 12(n - 2) + 8 = iV 3 . (125) 

It takes N 3 numbers to represent the Fourier transform of a cube N x N x N of real numbers, in spite 
of the fact that the Fourier transform is complex. Indeed, it is common for Fast Fourier Transform (FFT) 
subroutines to take a tridimensional array of number and to overwrite that array with the Fourier transform. 
This is the case for the FFT subroutines of Numerical Recipes (Press et al 1992). When these subroutines 
compute the Fourier transform of a cube N x N x N , the results are written in the same cube, and the actual 
numbers stored in the cube are precisely the 5 XXX 's, 8 XX 's, S x 's, and 5 U 's derived in this section. Hence it is 
possible to generate the density harmonics directly in fc-space, using the above formulae, store these numbers 
at the appropriate locations inside an N x N x N cubic array and then invoke the Numerical Recipes inverse 
FFT subroutines to generate the density field. 

The Numerical Recipes convention for storing the density harmonics is the following: Consider a 3D 
array A, with indicies running from to N — 1. We loop over all modes located in the first octant in fc-space: 
< l,m,n < N/2. 

1. For modes with < I, m, n < N/2 (the general case), the numbers 5 eee , d eeo , . . ., 5 000 are stored in a 
2x2x2 cube located at A(21,2m,2n), A(21,2m,2n+1), . . ., A(21+l , 2m+l , 2n+l) . 

2. For modes with I = 0, the numbers 6 ee , S eo , 5 oe , S QO are stored at A(0,2m,2n), A(0 , 2m, 2n+l) , 
A(0,2m+l,2n), and A(0 , 2m+l , 2n+l) , respectively. For modes with I = N/2, they are stored at 
A(l,2m,2n), A(l,2m,2n+1), A(l ,2m+l ,2n) , and A(l,2m+l,2n+l), respectively. This is easily gener- 
alized to the other faces (m = 0, m — N/2, n = 0, n = N/2). 

3. For modes with I — m = 0, the numbers S e and S are stored at A(0,0,2n) and A(0,0,2n+1), respec- 
tively. For modes with I = 0, m = N/2, they are stored are stored at A(0,l,2n) and A(0,l,2n+1), 
respectively. This is easily generalized to the other edges. 

4. The number 5 U is stored in a 2 x 2 x 2 cube located in the corner of the array, at A(0,0,0), A (0 , , 1 ) , 
. . ., A(l,l,l). Then, the value of A(0,0,0) is set to 0, since this represents the mode (0,0,0). 

4. CALCULATION OF THE DENSITY FIELD 

4.1. Eulerian Representation 

With the expressions derived in §2, we have all the ingredients necessary to compute the density field 
on a grid. The steps are the following: 

1. Choose a particular power spectrum P(k), a box size Lbox, and a grid size N. This determines the 
fundamental wavenumber fco = 2tt / 'Lbox- The allowed modes are given by equation (31). 

2. For all modes with l,m,n ^ 0, N/2, group these modes in groups of 8, and for each group, calculate 
the quantities S eee , S eeo , . . ., 5 000 using equations (70)-(77) and (80). The quantities R\, I\, . . ., I4 are 
determined from equations (60) and (61). 
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3. For all modes located on a face (one of the indicies l,m,n equal to or N/2), group these modes in 
groups of four, and for each group, calculate the quantities S ee , S eo , S oe , S 00 using equations (98)-(101) 
and (102)-(105). 

4. For all modes located on an edge (two of the indicies l,m,n equal to either or N/2), group these modes 
in groups of two, and for each group, calculate the quantities S e and S , using equations (120)-(101). 

5. For the modes located in a corner (all indicies I, m, n equal to either or N/2), calculate the quantity 
S u using equation (124). For the mode (I, m, n) = (0, 0, 0), replace the value of 6 U by 0. 

6. Store the quantities S eee , ■ ■ ■, S 000 , S ee , S eo , S oe , S 00 , S e , S , 5 U in a 3D, N x N x N array, at the proper 
locations. These depends on the convention used by the FFT subroutines. 

7. Compute the inverse FFT of the 3D array. The result will be the density field S(r). 

4.2. Lagrangian Representation 

In the Lagrangian representation, the density field is represented by a distribution of equal-mass par- 
ticles. We start by laying down N p x N p x N p particles on a cubic grid with grid spacing d = L^ ox /N p 
inside the computational volume Vb ox - The mass of the particles are given by M = pVbox/N 3 , such that the 
mean density inside the box is equal to the mean background density p. The particles are then displaced 
in order to represent the initial density field. This approach is valid only in the linear regime, defined by 
\S(r)\ <C 1. In this regime, the particle displacements Ar are significantly smallar than the initial separation 
d between the particles, and can be computed using the Zel'dovich approximation (Zel'dovich 1970). The 
displacements are given by 

k 

where Vj is the position of particle j before it is displaced, and rj + Arj is the position after. Notice that 
the reality condition (5) ensures that Ar^ is real. 

This method is straightforward, but can rapidly become unpractical. Equation (126) involves a sum 
over N 3 modes, and that sum must be performed for each of the N 3 particles. Since in typical cosmological 
simulations N p is chosen to be N/2, the number of operations scales like TV 6 . If it takes 5 minutes to set up 
initial conditions with 64 3 particles, it will take 5.3 hours for 128 3 particles, 2 weeks for 256 3 particles, and 
2.5 years for 512 3 particles! An alternative method, which scales like iV 3 , was proposed by Efstathiou et al. 
(1985). Essentially, this approach uses the fact that the displacement of each particle is proportional to its 
peculiar acceleration, that can be calculated with a N-body simulation algorithm such as PM (Particle-Mesh) 
or P 3 M (Particle-Particle/Particle-Mesh). With this method, the number of operations scales roughly as N 3 . 
I refer the reader to Hockney & Eastwood (1981) and Efstathiou et al. (1985) for details. It is worth noting 
that, unlike equation (126), the approach of Efstathiou et al. (1985) is approximative, and the high-frequency 
modes, near the Nyquist frequency k = (N p /2)k , are often poorly represented by the particle distribution. 

5. FILTERING 

Our next task in to filter the density field at some scale s. The choice of scale must obey two conditions: 
s» A and s <C L^ox- The first condition is required by the discreteness of the grid and the second by the 
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assumption of periodic boundary conditions. The filtered density field S s (r) at scale s is given by 



f 5{v')K s {v-v')d\' , (127) 

JVbox 



where if s is the filter function. We will use a Gaussian filter given by 

-x 2 /2s 2 



This filter function satisfied the normalization condition, 



/ K s {r)d 3 r = l, (129) 
Jv box 



as long as s <C £bc 



It is well known that filtering in real space is equivalent to a multiplication in fc-space. However, it is 
useful to redo the derivation, to ensure that we have all the correct factors of 2n, N 3 , and so on. First, we 
express the filter as an inverse Fourier transform, 



K s (r -r') = i^ KMe-*' <*-'"> . (130) 

k' 

We substitute equations (33) and (130) in equation (127), and get 

w = 4e / E^ k ) e ~* k ' r E^( k >~ 4k '' (r ~ r ' w 

= 4EE^ k )^( k ')^ lkr / e t(k '- k)r 'dv. (isi) 

k k' JV b°x 

(132) 

The integral is equal to Vb x<5k,k' (see eq. [12]), and we use the Kronecker S to eliminate the sum over k'. 
We get 

^( r )-^E^ k )^( k ) e " k ' r - ( 133 ) 

k 

We now need an expression for K s (k). This function is the Fourier transform of the filter, 

k s (k) = *.(*)e*- = 7^7^ E e ^ 2/2s2 e Jk ' x ■ (134) 

We rewrite this expression as 

N * \ -x 2 /2 s \ ik-x f V bc 



- 22 - 



The factor Vbox/-Y 3 = A 3 represents the volume element around each point x in the N x N x N grid. Since 
we assume s » A, we can approximate the sum as an integral over the volume of the box (or, equivalently, 
regard the sum as a numerical approximation for the integral). Hence, 

Since we assume periodic boundary conditions, we are free to locate the origin anywhere inside the box. For 
instance, we can locate it in the center of the box. Since s <C Lbox, the integrant becomes negligible at the 
edge of the box. We can then extend the integration domain to all space, 

where this expression no longer assumes boundary conditions. The integral in equation (137) can be found 
in any textbook of Fourier transforms, 

f e - X Z/2s* ei k- Xd 3 x ^ (27r) 3/2 s 3 e -(fcs) 2 /2 . (138) 

Jail space 

Hence, 

^(k) = #^- (fcs)2/2 . (139) 

Vbox 

We substitute this expression in equation (133), and get 

^ r )^E^ N2/v,k ' r - ( 14 °) 

k 

Hence, to obtain a filtered density field, we generate the density harmonics using the method described in 
§3, and then multiply them by the factor e~ k s / 2 before taking the inverse Fourier transform. 



6. SUMMARY 

This paper presents in great detail the techniques used for generating Gaussian density fields. These 
techniques are well-known among experts in cosmological numerical simulations, but the specific details of 
the implementation are often difficult to find in the literature. Also, the notation tends to vary significantly 
from one author to another. The consequences is that any new researcher moving into this field has to cither 
spend a great deal of effort rederiving all the technical details, or else rely on existing codes and use them as 
black boxes. The goal of this document is to improve the situation by presenting in a comprehensive form 
the basic theory behind the generation of Gaussian random fields. 
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A. FOURIER TRANSFORM OF THE TOP HAT 

Consider the following integral, 



i= f dV~ 4ky , 

Jsph(O) 



(Al) 



where the domain of integration is a sphere of radius R centered at the origin. We consider a spherical 
coordinate system centered at the origin, with the z-axis pointing in the direction of k. Equation (Al) 
becomes 

r 2ir /-7T rR 

^—iky cos 9 



/* Z7T PIT Pit 

1=1 dej) / d9 I dy{y 2 sirie)e- 
Jo Jo Jo 



(A2) 



where k = |k|, y = |y|, and 9 is the angle between k and y. The integrations over <f> and 9 are trivial. We 
get 

-iky COS0-I 71 " f-R _ [ piky _ p -iky "| fR y sill ky 



I = 2tt[ dyy 2 
Jo 



iky 



= 2ir 
o Jo 



[ dyy 2 
Jo 



iky 



The integral over y is now trivial. We get 



47r 4ttR 3 
I = — (sin kR — kR cos kR) = — — (sin u — u cos u) , 
k A u A 



(A3) 



(A4) 



where u = kR. 
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